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ABSTRACT  XI 

A  study  was  conducted  to  numerically  treat  a 
mixture  composed  of  a  gas  and  solid  particles  in  a  supersonic* 
axi symmetric  nozzle.  The  governing  equations  are  a  set  of  eight 
first  order,  quasi-linear ,  partial  differential  equations. 

Seven  of  these  equations  are  of  the  hyperbolic  type  when  the 
flow  is  supersonic  (based  on  the  frozen  speed  of  sound  in  the 
gas)  and  can  be  solved  by  the  method  of  characteristics.  The 
eighth  equation  (the  particle  continuity  equation)  is  rewritten 
as  an  integral  equation  to  be  solved.  The  resulting  seven 
compatibility  equations  and  the  seven  characteristic  equations 
(only  four  are  distinctj  the  t.  gas  Mach  lines  and  the  gas 
and  particle  streamlines)  are  solved  by  the  modified  Euler 
predictor-  corrector  algorithm.  These  equations  were  programmed 
for  an  IBM  1130  computer.  A  sample  nozzle  calculation  is  given 
and  compared  with  the  one-dimensional  calculations.  These 
results  indicate  that  the  program  is  working  correctly.^ 


I.  INTRODUCTION 

The  purpose  of  the  present  investigation  was  to 
numerically  treat  a  mixture  of  a  gas  and  solid  particles  in 
an  axisymmetric ,  supersonic  nozzle.  The  study  vas  primarily 
intended  for  industrial  particle  micronization  processes.  In 
this  process  solid  particles  are  fluidized  and  entrained 
by  a  high  pressure  air  flow.  The  mixture  is  then  expanded 
through  a  converging-diverging  nozzle.  The  solid  particles  are 
accelerated  during  this  expansion  and  impact  a  target  downstream 
of  the  nozzle  exit  plane.  If  the  particles  have  sufficient 
momentum  they  are  broken  upon  impact.  This  process  is  repeated 
until  the  desired  particle  sizes  are  obtained. 

The  interest  in  the  present  study  was  not  so  much 
to  determine  two-dimensional  particle  velocities  but  rather  to 
be  able  to  investigate  non-uniform  particle  distributions  in 
each  nozzle  plane.  One  of  the  inherent  restrictions  in  one¬ 
dimensional  gas  solid  -  particle  analyses  is  that  the  particles 
are  uniformly  distributed  in  each  nozzle  plane.  However, 
experimental  studies  definitely  show  that  most  of  the  particles 
are  concentrated  near  the  nozzle  center  line.  To  treat  the 
problem  of  nonuniform  particle  distribution  it  is  necessary 
to  formulate  the  problem  in  two  dimensions. 


Riethmuller1  h-vs  done  an  extensive  literature  survey 
on  ga3-solid  particle  flow  analyses  and  experiments.  Therefore, 
only  articles  that  are  directly  pertinent  to  this  analysis 
are  referenced  in  this  report. 


XI.  GOVERHlNG  EQUATIONS 

The  supersonic  motion  of  most  compressible 
fluids  encountered  in  nozzle  expansion  flow  can  be  accurately 
described  by  the  governing  equations  of  an  invieciu  fluid. 

The  basic  assumptions  which  constitute  such  a  gas  dynamic  model 
are  :  (1)  continuum,  (2)  steady,  (3)  inviscid,  (k)  adiabatic, 
and  (5)  the  gas  composition  is  frozen.  The  resulting  partial 
differential  equations  can  be  treated  by  the  veil  known 
technique  of  the  method  of  characteristics.  For  axisymnetric 
or  two-dimensional  nozzles  this  method  transforms  the  partial 
differential  equations  into  a  system  of  differential  equations 
which  are  valid  along  certain  characteristic  directions  in  the 
flow  field. 


The  transformed  equations  are  much  simpler  to  solve 
and  this  fact  led  Kliegel2  to  attempt  a  similar  analysis  for 
a  mixture  of  a  gas  and  solid  particles.  The  major  assumption 
necessary  f-'r  such  a  gas-solid  particle  analysis  is  that  the 
particles  hehave  as  a  continuum.  Clearly,  this  assumption  is 
physically  unrealistic.  That  it  yields  good  engineering  results 
is  another  question.  Kliegel  has  shown  that  the  results  are 
in  good  agreement  with  observation  for  rocket  engine  applications 
for  moderate  loading  ratios  and  small  particle  sizes.  Similar 
comparisons  will  he  required  for  particle  micronization  processes 
to  determine  the  applicability  of  this  assumption.  However, 
it  would  appear  that, even  for  particle  micronization  processes 

this  assumption  should  give  useful  engineering  results  for 
moderate  loading  ratios  of  small  particles. 

The  equations  governing  the  steady  axisymmetsric  flow 
of  a  mixture  of  a  gas  and  solid  particles  are  derived  by 
Hoffman  and  Thompson3.  The  assumptions  necessary  for  these 
derivations  are  :  (1)  the  gas  is  inviscid  except  for  its 
interactions  with  the  solid  particles,  (2)  the  mixture  mass  and 
mixture  energy  of  the  system  are  constant,  (3)  the  volume 


3 
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occupied  by  the  solid  particles  is  negligible,  (4)  the 
therisal  notion  of  the  solid  particles  ia  negligible,  (5) 
the  solid  particles  do  not  interact,  (6)  the  drag  and  heat 
transfer  characteristics  of  an  actual  particle  shape  and  the 
size  distribution  of  particles  can  be  represented  by  spherical 
particles  of  a  single  size,  (7)  the  internal  temperature  of 
each  solid  particle  is  uniform,  (8)  energy  exchange  between 
the  gas  and  the  solid  particles  occurs  by  convection  only, 

(9)  the  only  force  acting  on  the  solid  particles  is  the 
viscous  drag  fo'-es,  (10)  no  mass  transfer  between  the  gas 
and  the  solid  particles,  and  (11)  no  phase  change. 


Based  on  these  assumptions,  the  following  equations 
govern  the  gas-solid  particle  flow  : 
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The  drag  coefficient  and  Nusselt  nuaber  for  a 
particle  are  calculated  from  equations  given  by  Neilaon 
and  Gilchrist7 


X*  ^  •  ; 

— -r  ^ 


•  ■  O 

t-  T-“ 

r  - 


28  Re”  °*85  +  0 . 48 


N  «  2.0  +  0.6  Re°*s  Pr G  ‘  3  3  3 
u 


(11-12) 

(11-13) 


vhere  the  Reynolds  nuaber  for  c,  particle  is  defined 
as  follows  : 

pD(W  -  S  ) 
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(II-U) 


III.  APPLICATION  OP  THE  METHOD  OF  CHARACTERISTICS 


In  this  section,  the  techniques  of  the  method  of 
characteristics  vill  be  employed  to  obtain  the  characteristic 
and  compatibility  equations  for  the  flow  fi  ,ld  variables. 

The  flow  field  governing  equations.  Equations  ( XI— 1 ) 
through  (II-6)  can  be  written  in  the  following  general  form  : 

L.  *  a..z.  +  b..z.  +  c.  «  0  (i»j  *  8)  (III-1) 

i  ij  jx  ij  jy  i 

where  z  represents  the  eight  dependent  variables  u,v,p,p,u^, 

v  ,  p  and  T  and  the  x  and  y  subscripts  denote  partial 
P  P  P 

differentiation.  These  eight  equations  can  be  combined  to 
form  a  single  differential  operator  by  employing  arbitrary 
multipliers  and  summing.  Thus 

L  *  o.  L.  »  0  ( i  «*  1,...,  8)  (III-2) 

i  x 

where  are  the  arbitrary  multipliers.  In  this  section,  the 
convention  of  indicati' 3  a  summation  be  repeated  indices  will 
be  used.  The  partial  differential  equation.  Equation  {111-2} , 
can  be  rewritten  in  the  form  of  an  ordinary  differential 
equation  under  certain  conditions.  Thus 

(a-.o.)dzj  +  c.o.dx  =  0  (i,j  »  1,...,  8)  (III-3) 

ij  x  1  1 

if  and  only  if  the  following  equations  are  valid  : 

a.  .0.  dy/dx  *  b..<j.  (i»<j  *  ^  *  ®)  vill— 4} 

1  U  1 

Equations  (III-4)  are  eight  independent  equations  for  dy/dx. 

Equations  (JII-4)  are  used  to  determine  the  unknown  multipliers 

a.  which  can  then  be  substituted  into  Equation  (ni~3)  to  yield 
1 

[  the  compatibility  equations.  Rearranging  Equations  (111-4)  to 

|  solve  for  yields  the  following  equations  : 
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If  this  system  of  equations.  Equations  (III-5),  ia  to  have 
a  solution  other  than  the  trivial,  i.e.,  all  *  0,  then 
the  determinant  of  the  coefficients  of  must  be  aero.  Thus, 


[  a. .  dy/dx  -  b .  ,  ]  3  0 
L  10  lj  J 


( i  ,  j  «  1.....8)  ( III-6) 


Nov,  dy/dx  can  be  obtain sd  from  Equation  (III-6).  The 
resulting  dy/dx  can  then  be  substituted  into  Equations  flll-U) 
to  determine  relationships  in  terms  of  the  multipliers  o^. 

With  Oj,  known,  the  final  form  of  the  compatibility  equations. 
Equations  (III-3),  is  obtained. 

The  governing  equations  for  the  flow  field  are 
repeated  below  for  convenience. 
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The  determinant,  Equation  (III-6)  can  now  be  written  as 
follows  : 


pdy/dx  -  p  0 


pSj  0  dy/dx  0 


d  dv/dbac  -  p  0 
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PpS2 


V* 


0  p  cSo  0  I 
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(111-15) 


where  S.  -  (udy/dx  -  v)  and  S2  *  (u  dy/dx  -  v  }. 

P  P 

Since  the  upper  right  quarter  of  the  above  determinant  is  filled 
with  zeros,  the  expansion  of  the  determinant  reduces  to 


[  G  1 


[  f  ]  *  ° 


( III-16) 


where  G  represents  the  upper  left  quarter  of  the  determinant 
and  P  represents  the  lower  right  quarter  of  the  determinant 
as  shown  below  : 
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Eq.uftti.or.  (III-16)  c*n  b*  zero  by  either  G  or  P  being  zero. 
Setting  these  two  determinants  to  zero  yield*  the  following 
expressions  : 

[0]  -  ( udy/dx  -  v>2  i Uy/a*>2  I  «!  •  *!  I  -  2  «ayM* 

♦  (Ta  -  „2)]  .  o  (HI-17) 


ftnd 


[pl  *  (updy/dx  -  vp)4  *  0 

Therefore,  Equation  (HI-1 6)  becomes  : 

(u  dy/dx  -  v  )4  (udy/dx  -  v)2  [  (dy/dx)2[  u2  -  a2] 
P  P 


(III-18) 


2  uvdy/dx 


+  (v2  -  a2)  j  -  0  (III-19) 


The  characteristic  curves  are  found  by  solving  for  dy/dx 
which,  for  this  system  of  equations,  are 


dy/dx  a  v/u 

uv  i  a2  /  m'2'~  1 


dy/dx 


u2  -  a2 


dy/dx  *  vp/up 


(III-20) 
(III-21 ) 
(III-22) 


!»  teras  of  t*e  Mac*  angle  a  .  t*.  8«  6  th« 


LXi  verwa  - - «  .  .  « 

particle  flow  angle  8p,  Equations  (111-20),  (HI-21 )  an 
(III-22 )  can  be  rewritten,  respectively,  as 


dy/dx  *  tan  e 


(III-23) 


dy/dx  *  tan  (6  ±  o) 


(III-2h) 
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dy/dx  *  tan  6p 
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Thus,  the  characteristic  curves  sre  the  gas  streamlines 
appearing  twice,  the  Mach  lines  each  appearing  once  and  the 
particle  streamlines  appearing  four  times.  Equations  (III— 21  ) 
are  real  only  when  M  >  1,  vhereas.  Equations  (III-2Q)  and 
(III-22)  are  always  real.  Therefore,  the  system  of  governing 
equations.  Equations  ( I I I — T )  through  (III-lU)  is  hyperbolic 
when  the  flow  field  is  supersonic.  Of  these  eight  characte¬ 
ristic  curves,  four  are  distinct;  the  two  Mach  lines,  the  gas 
streamlines  and  the  particle  streamlines. 


In  order  to  determine  the  compatibility  equations, 
the  unknown  multiplier  must  be  determined.  These  multipliers 
are  determined  by  simultaneously  solving  Equations  (Ill-k), 
using  Equations  (III-2Q)  through  (III-22)  for  the  slope  of 
each  characteristic  curve. 


Substituting  the  coefficients  a^  and  b.^  from 
Equations  (IXI-7)  through  (III-lU)  into  Equations  ? Ill— 1* ) 
yields  the  following  results  : 


a ldy/dx  +  (udy/dx  -  v)o2  =  0 


(III-26) 


-0}  +  (udy/dx  -  v)o3  3  0 


(HI-27) 


r2dy/dx  -  o3  +  (udy/dx  -  v)o4  3  0 


(111-28) 


(udy/dx  -  v)  (oj  -  a20t,)  3  0 


(III-29) 


o5dy/dx  ♦  (u^dy/dx  -  v?)o6 


(III-30) 


-oc  +  (u  dy/dx  -  v  )<j7  *  0 
3  P  P 


(u^dy/dx  »■  vp)o8 


( III-31  ) 


(III-32) 


(u  dy/dx  -  v  )o5  *  0 
P  P 


(III-33) 
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The  fora  of  the  general  compatibility  equation. 
Equation  (II 1—3 )  is  : 

[  pOj+  puo2  ]  du  +  [PUO3  ]  dv  +  [o2  *  u04  ]  dp 

+  [  uOi  -  a2ua4  ]  dp  +  £ppo5  *  ppUp°6^dup 

+  [  p  u  a7  ]  dv  +  [  cp  u  Co  ]  dT 
LPP7J  P  1  P  P  8  J  P 

+  [  up®5  J dpp  +  [  pvoj/y  +  Ap^(u  “  up^°2 


+Ap  (t  -  v  )aj  -  ABp  a4  +  p  v  a5/y  -  Ap  (u 
P  p  P  P  P  P 

-  Ap  (v  -  v  )a 7  +  -Ip  AC ( T  -  T)o8  ]  dx  *  0 
P  P  3  p  p  1 


»p)o6 


( III-3U ) 


For  the  gas  streamline  characteristic  curve, 

(udy/dx  -  v)  *  0.  Thus,  by  substituting  thi3  relationship  into 
Equations  ("11-26)  through  (III-33),  the  following  results 
are  obtained  : 

Oi  *  0 


a 2  *  U03/V 


( III-35 ) 


Or,  »  <J5  *  07  »  Og  *  0 

Substituting  Equation  (111-35)  into  Equation  (lII-3*»)  and 
regrouping  tera3  into  coefficients  of  the  two  arbitrary 
multipliers ,  03  and  o4,  yields  the  following  result  : 

[  pudu  +  pvdv  +  dp  +  Ap  (u  -  u  )dx  +  Ap  (v  -  v  )dy  ]  ^o3 

*■  Tl  TJ  u  ¥ 


+  [  udp  -  a2udp  -  ABp^dx  ]o4  *  0 


(111-36) 


Since  the  multipliers  o3  and  Oj,  are  arbitrary,  their  coefficients 
must  equal  xero.  Equating  the  coefficient#  of  these  multipliers 


-  1 1  - 


to  zero  yields  the  following  compatibility  equations  which 


ar*  valid  along  gas  streamlines  : 


pudu  +  pvdv  +  dp  +  Ap  (u  -  u  )dx  +  Ap  (v  -  v  )dy  »  0  (III-37) 

P  P  P  P 


udp  -  a2udp  -  ABp  d:t  *  0 


(111-38) 


Votings  that 


udu  +  vdv  *  WdW 


Equation  (III-37)  can  be  rewritten  as 


pS?dV  +  dp  +  Ap  [  (u  -  u  )dx-  *  (v  -  v  )dy  ]  »  0  (III-39) 

P  P  P 


Equations  (111-38)  and  (III-39)  are  valid  only  aldng  gas  stream¬ 
lines,  i.e.,  along  curves  defined  by 


dy/dx  *  tan  0 


( 111-1*0) 


Along  the  Mach  lines,  the  slope  dy/dx  is  determined 
from  either  Equations  (III-21)  or  (III-2U)  and  (udy/dx  -  v)  +  0. 
Therefore,  by  combining  Equations  (111-26)  through  (III-33) 
with  Equation  (III-31),  the  relationships  for  the  multipliers 


a j  for  determining  the  compatibility  equations  which  are  valid 


along  the  Mach  lines  are  obtained  and  are  as  follows  : 


a2Oi* 


-  a  2 


a2oi*  dy/dx 


(udy/dx  -  v) 
f  (ud> /dx  -  v) 2 


-  a2 (dy/dx) 2 


(udy/dx  - 


®5  “ 


(Ill-Ul ) 
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Since  no  relationship  is  obtained  for  c4„  thi >j  multiplier 
is  arbitrary.  The  compatibility  equations  ralid  along  the 
Mach  lines  are  obtained  by  substituting  Equation*  (111-1*1) 
into  Equation  (IXI-3^),  the  general  compatibility  equation, 
and  equating  the  coefficient  of  o4  to  tero,  which  yields 

( “  vdu  ♦  udr)  ±  1-  dp  +  — (udy  -  rdx) 


-  A(p  / p )  (u  -  u  )dy  +  A(p  /p)(v  -  ▼  )dx 
P  P  P  p 


~  AB( pp/pa2 ) ( udy  -  vdx)  »  0 


(111-1*2) 

Equation  ( 111-1*2}  can  be  rewritten  in  terns  of  8  and  a  as 


•  /i 

d8  ±(cota/oW2)dp  ±  £  -S-- ^-6' - (1  +  B/a2)}  * 

('c"ors(e"V  aj)dx  +  A(pp /pw2)  r  uPdy  "  V*  ^  “  ° 


( 111-1*3 ) 


where  the  ga3  Mach  lines  are  given  by 


dy/dx  *  tan  (8  ±  a) 


(111-1*1*) 


In  Equations  (111-1*3)  and  (111-1*1*),  the  upper  signs  refer  to 
left-running  M»ch  liner  and  the  lower  signs  refer  to  right¬ 
running  Mach  lines. 

For  the  particle  streamline  characteristic  cur re, 
(uyiy/dx  -  ▼  )  ■  0,  Thus,  by  substituting  this  relationship 
into  Equations  (III-26)  through  (III-33),  the  following  results 
are  obtained  : 


•  °2  *  °3  *  °4  ®  °5  “  0 


( 111-1*5) 
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Substituting  Equation  (111-45)  into  Equation  (111-34)  and 
regrouping  terms  into  coefficients  of  the  three  arbitrary 
multipliers,  og,  07  and  o8  yields  the  following  result  : 


f  u  du  -  A(u  -  u  )dx  1 0 c  ♦  [  u  dv  -  A(v  -  v  )dx  1 07 

k  p  P  P  J6  LPP  PJ 


+  f  ca  dT  f  AC ( Tp  -  T)dx  ]aB 

p  p  3  -*  v 


(111-46) 


Again,  the  multipliers  og ,  07  and  o8  are  arbitrary.  Therefore, 
their  coef ficients  must  be  zero.  Equating  the  coefficients  of 
these  three  multipliers  to  zero  yields  the  following  compati¬ 
bility  equations  which  are  valid  along  the  particle  streamlines 


u  du  -  A(u  -  u  )dx  «  0 

p  p  p 


(111-47) 


u  4v  -  A(v  -  v  )dx  *  0 
P  P  P 


( III-48 ) 


cu  dT  +  ~  AC ( T  -  T)dx  *  0 
p  p  3  P 


(111-49) 


where  the  particle  streamline  is  given  by 


dy/dx  =  v  /u 
P  P 


(III-50) 


It  should  be  pointed  out  that  only  seven  compatibility 
equations  are  obtained  in  this  analysis,  i.e.,  no  equation  is 


obtained  for  solving  for  the  particle  density,  This  can 


be  explained  by  noting  that  05  is  zero  for  all  characteristic 
curves.  Therefore,  the  fifth  governing  equation.  Equation  (III-ll), 
the  particle  continuity  equation  does  not  inter  the  analysis 
as  treated  by  the  method  of  characteristics.  If  Equation  (III-ll) 
is  rersoved  from  the  set  of  governing  equations,  seven  charac¬ 
teristic  curves  and  seven  compatibility  equations  are  obtained, 
as  should  be.  In  order  to  have  a  complete  set  of  equations  it 
is  necessary  to  include  Equation  (III— 11)  or  some  equivalent 
equation  in  order  to  have  a  relationship  for  p^.  In  the  present 
analysis  an  integral  equivalent  equation  for  particle  continuity 


*- — 


V 


is  employed.  The  equation  used  is  as  follows  : 


fy 

m  »  2*  I  yp  (u  dy/dx  -  v  )dx  (IT  1—5 1 ) 

P  P  P  P 

v 

Equation  ( III— 51 )  is  used  to  calculate  an  average  particle 
density  at  each  point  vn  the  flow  field. 

It  should  be  reemphasized  that  treating  the  particles 
as  a  continuum  is  an  approximation  in  the  present  two- 
dimensional  analysis  just  as  in  previous  one-dimensional 
analyses . 


In  summary,  the  governing  partial  differential 
equations  have  been  transformed  to  a  system  of  differential 
equations  valid  only  along  certain  characteristic  curves. 
These  are  : 

Gas  Streamlines 


dy/dx  •>  tan  8 

pWdW  +  dp  +  Ap  [  (u-u)dx+(v-v  )dy] 

P  p  P 

udp  -  a2udp  -  ABp  dx  =  0 

P 

Gas  Mach  Lines 
dy/dx  -  tan  (6  ±  a) 

d8  ±(cota/pW2)dp  ±  [  (1  +  B/a2)] 

+  A(pp/pW2)  [  u^dy  -  v^dx  ]  *  0 


( 111*52) 
*  0  (III-S3) 

(III-SU) 


(111*55) 

- Sf;V-T 

cos  (8  ±  a) 

( 111-56) 


Particle  8tre» val ines 


dy/dx  *  tan  8p 

u  du  -  A(u  -  u  )dx  *  0 
P  P  P 

u  dv  -  A(v  -  v  )dx  *  0 
P  P  P 

u  cdT  +  §  AC ( T  -  T)dx  «  0 
P  P  3  p 

and  the  particle  continuity  equation 

ry 

m  *  2u  yp  (u  dy/dx  -  v  )dx 


(111-57/ 

(111-58) 

(HI-59) 

(III-60) 


(III-61 ) 
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IV.  NUMERICAL  SOLUTION  TECHNIQUE 

The  formulation  of  the  problem  of  treating  a 
supersonic  flow  field  consisting  of  &  mixture  of  a  gas  ana 
solid  particles  represents  only  one  part  of  the  present 
investigation.  In  order  to  make  use  of  the  equations  derived 
in  the  previous  section,  it  is  necessary  to  translate  the 
equations  of  motion  of  Section  III  into  a  numerical  algorithm. 
This  section  outlines  the  solution  procedure  which  was 
programmed  for  the  IBM  1130  computer  at  the  von  Karman  Institute. 

A3  was  discussed  in  the  previous  section,  the  gas 
Mach  lines  (two  of  the  characteristic  curves)  are  real  only 
when  the  gas  Mach  number  (ratio  of  the  gas  speed  and  the  frozen 
speed  of  sound)  is  greater  than  one.  Therefore,  the  present 
scheme  is  only  applicable  in  regions  of  the  flow  which  are 
•■’upersonic.  This  means  that  in  order  to  calculate  the  flow  in 
a  nozzle  by  the  method  of  characteristics  one  must  know  the 
flow  properties  in  the  region  downstream  of,  but  near,  the  gas 
sonic  line.  In  order  to  obtain  the  necessary  initial  conditions 
for  a  gas-solid  particle  flow  the  complete  subsonic  and  transonic 
flow  field  must  be  solved.  Since  thin  problem  has  not  been 
solved  to  date,  it  is  necessary  to  use  approximate  methods 
to  obtain  the  gas  and  solid  particle  properties  from  which  to 
start  the  solution  procedure.  It  is  clearly  evident  that  much 
more  work  must  be  done  in  the  subsonic  and  transonic  flow 
regions  before  the  present  supersonic  analysis  can  be  effecti¬ 
vely  utilized,  i.e.,  to  get  realistic  results  one  must  specify 
realistic  initial  conditions. 

The  numerical  algorithm  does  not  depend  on  the 
initial  conditions.  Therefore,  as  better  methods  for  predicting 
the  initial  conditions  become  available,  they  can  be  input 
into  the  present  analysis  directly  as  input  data. 

The  numerical  algorithm  used  is  a  modified  Euler, 
predictor-corrector  technique  where  averaged  coefficients  are 
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i. 

j, 

V 

used.  The  transformed  equations.  Equations  (II.I-52)  through  \ 

(III-61)  are  forward  differenced.  To  illustrate  the  technique,  i 

Figure  1  presents  an  interior  mesh  point  gride  The  four  distinct  . 

characteristic  curves  are  shown  with  the  point  nustbering  scheae  j 

as  used  in  the  computer  program.  Point  3  is  the  point  where 
the  solution  is  sought.  Points  1  and  2  are  previously  obtained, 
solution  points  or  given  points  on  the  initial  value  line. 

The  prop  ies  at  points  4  and  5  are  known  (by  interpolation) 
once  the  vocation  of  each  point  is  determined.  The  solution 
procedure  is  as  follows  : 

(1 )  Predict  the  x  and  y  direction  of  Point  3  by  solving 
the  two  equations.  Equations  (III-55)*  In  finite  difference 
form  this  yields  : 

*3  *  ^2  '  n  -  *2&?.3  +  XlSi3)/\Sl3  -  h23) 
y 3  ■- yi  *  s13(x3  -  xj) 


where 


H2  3  ~  ^ 

[  tan 

(  8  2  +  '*2  ^ 

+  tan 

(e3 

+  a  3 )  J 

1 

'-13  *  '2 

[  tan 

(01  "  Qj ) 

*■  tan 

( 9  3 

+  a3)  j 

it  should  be  noted  that  the  second  term  within  the  brackets 
in  the  above  two  equations  are  set  equal  to  the  first  term 
for  the  predictor. 

(2)  Solvs  for  the  gas  flow  ang3e  and  the  pressure  by  simulta¬ 
neously  solving  the  two  equations,  Equations  (111-56). 

(3)  Calculate  the  position  of  Points  h  and  5  by  the  particle 
streamline  equation.  Equation  (111-57)  and  the  ge.s  streamline 
equation.  Equation  (111-52),  respectively. 
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(U)  Interpolate  for  the  properties  at  points  U  and  5* 

(5)  Calculate  the  gas  velocity  and  the  gas  density  at  P»int  3 
from  Equations  (III-53)  ssd  (IIX-5M. 

(6)  Calculate  the  particle  velocity,  particle  flow  angle, 
particl®  temperature,  and  particle  density  from  Equations 
(111-53)  through  (III-61). 

(7)  Nov  that  all  properties  are  known  (predicted  values), 
evaluate  all  equation  coefficients,  average  the  coefficients, 
and  repeat  steps  (i)  through  (6). 

This  algorithm  is  second  order  accurate.  To 
iterate  again  accomplishes  nothing.  If  the  accuracy  is  not 
sufficient  change  the  grid  spacing-by  increasing  the  number 
of  points  on  the  initial  value  line- and  redo  the  calculations. 

The  same  technique  is  used  on  the  vail  boundary 
except  that  on  the  vail  the  solution  point  is  obtained  by  the 
intersection  of  the  left  characteristic  and  the  vail  crntour. 
Since  the  gas  flow  angle  is  the  vail  angle,  the  flov  angle  is 
known.  Thi3  reduces  the  unknowns  by  one,  which  is  necessary 
since  the  right  characteristic  compatibility  equation  cannot 
be  used. 

Tbe  same  procedure  exists  at  the  centerline. 

Here  the  gas  flow  angle  is  zero  (known,  1  priori),  and  the 
left  characteristic  compatibility  equation  is  not  used. 

To  check  the  program  out,  a  sample  ease  was 
executed  for  which  previous  one-dimensional  results  were 
available.  The  initial  value  li-.c  properties  were  taken 
from  the  one-dimensional  results.  The  resulte  were  then 
compared  at  the  nozzle  exit,  9  cm.  downstream  of  the  initial 
value  line.  The  following  table  presents  the  results  a r. 
obtained  from  the  two  programs.  The  loading  ratio  was  1.5, 


i^e.,  the  particle  flow  rate  was  1.5  times  the  gas  flow  rate. 
Particle  diameter  was  500  microns. 


Gas  Vel&eity 
Particle  Velocity 
Ga3  Temperature 
Particle  Temperature 
Gas  Hach  number 
Pressure 

The  run  time  for  the  two-dimensional  calculation 
on  the  IBM  1130  was  approximately  9  minutes. 

The  above  comparison  indicates  that  the  program  is 
working  correctly  and  that  it  can  be  used  effectively  for 
perticle-gas  problems. 


Present  Analysis 

2  -  d  1  -  d 

t6?.0  m/sec  U67 .6  m/sec 

213.2  m/sec  21h.7  m/sec 

159-7  °K  159. 7°K 

295.7  °K  291 .9  °K 

1.950  1.810 

1.86x106n/m2  1.88xl06  a/m 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  results  of  the  comparison  of  the  present 
two-dimensional,  supersonic  analysis  with  the  one-dimensional 
analysis  indicate  that  the  program  is  working  correctly 
and  that  it  can  be  a  useful  tool  for  particle  micronization 
processes.  The  present  analysis  has  the  advantage  that 
non-uniform  particle  distributions  can  be  treated  if  the 
necessary  initial  data  are  available.  The  program  has  been 
written  in  modular  form  so  that  any  required  or  desired 
program  changes  are  relatively  simple. 

Further  studies  are  needed  in  two  areas.  These 
are  :  ( 1 )  a  working  subsonic  and  transonic  analysis  to  take 
advantage  of  the  potential  applications  of  this  analysis 
and  (2)  much  more  effort  is  required  to  obtain  governing 
equations  that  are  more  physically  realistic.  The  use  of  the 
particle  continuum  assumption,  obviously,  has  restrictions. 
Experimental  measurements  are  necessary  to  determine  the 
region  where  this  assumption  breaks  down.  It  should  be 
emphasized  that  one-dimensional  analyses  suffer  from  the  same 
problem  and,  therefore,  if  for  particle  micronization  processes 
the  continuum  assumption  does  not  provide  engineering 
results  neither  approach  can  be  used  to  predict  particle 
behaviour  in  the  nozzle. 
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APPENDIX  A 
PROGRAM  DESCRIPTION 

The  computer  program  was  written  in  the  FORTRAN  XV 
language  for  the  IBM  1130  computer  at  the  von  Karaan  Institute. 

The  program  consists  of  the  MAIN  program  and  four  subroutines. 
These  subroutines  are  designated FL0WP,  M0CP,  M0CB  and  C0EFP.  The 
structure  of  the  program  is  illustrated  in  Figure  A-1.  A  brief 
description  of  the  program  is  as  follows  : 

MAIN  -  This  program  call3  the  flow  field  logic  program 
or  GALLS  EXIT  as  determined  by  whether  SWITCH  10  oa  the  console 
is  on  or  oix  • 

FL0WP  -  This  subroutine  controls  the  logic  in  calcu¬ 
lating  the  flow  field  in  an  axisymmetric  nozzle.  All  input  data 
and  all  output  options  are  accomplished  by  this  subroutine. 
Subroutine  FL0WP  calls  M0CP  for  the  solution  at  an  interior  mesh 
point  or  M0CB  for  the  solution  at  the  wall  boundary  or  the  center¬ 
line  boundary.  Two  arrays  store  the  solution  along  a  bach  left 
characteristic  and  the  left  characteristic  being  solved. 

M0CP  -  This  subroutine  uses  the  modified  Euler  predic¬ 
tor  corrector  technique  with  averaged  coefficients  to  solve  for 
an  interior  nesh  point.  Subroutine  C0EFP  is  called  to  evaluate 
all  the  necessary  coefficients  used  in  the  finite  difference 
equations . 

M0CB  -  This  subroutine  is  nearly  identical  to  the 
above  subroutine  except  the  solution  point  is  on  the  boundary. 

C0EFP  -  This  subroutine  evaluates  the  coefficients 
of  the  finite  difference  equations  and  is  called  from  both 
M0c;?  and  M0CB. 
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APPENDIX  B 


PROGRAM  INPUT 

The  input  data  for  the  program  has  been  kept  to  a 
minimum  to  keep  the  program  simple.  The  following  data  are 
required  for  program  execution  : 

CARD  1  (FORMAT  3F10.0) 

I- 10  GAMMA,  ratio  of  the  specific  heats  for  the  gas  (air,  y*1.U) 

II- 20  RGAS ,  the  gas  constant  (air,  R=29.27) 

21-30  PR,  the  gas  Prandtl  number  (air,  Pr*0.7) 

CARD  2  (FORMAT  4F10.0) 

I- 10  DIAM,  particle  diameter  (meters) 

II- 20  RHOP,  particle  density  (kg/m3) 

21-30  CPART ,  particle  specific  heat 

3 1-4o  PFLRT ,  particle  flow  rate  (kg/sec) 

CARD  3  (FORMAT  12 ,8X . 12 5CX , 12 ,8X , I 1 ) 

I- 2  NPTSL,  the  number  of  points  for  which  data  is  specified 

on  the  initial  value  line  (MAX*5) 

II- 12  IMAX,  the  maximum  number  of  left  characteristics  to  be 

calculated 

21-22  IOUT ,  the  number  of  the  left  characteristic  where  print 
out  begins.  If  all  characteristics  are  desired  I0UT=0. 

31  NU,  not  used 

The  following  two  cards  are  repeated  for  each  of  the 
NPTSL  points.  Two  cards  are  required  per  point.  The  first  point 
input  is  the  wall  point,  the  last  point  is  the  centerline  point. 

CARD  4  (FORMAT  ?F10.0) 

I- 10  x-location  of  the  point  (meters) 

II- 20  y-location  of  the  point  (meters) 

2Y-30  Gas  Mach  number 
31-40  Gas  velocity  (m/sec) 

41-50  Gas  flow  angle  (degrees) 

51-60  pressure  (n/m2) 

61-70  gas  density  (kg/mJ  e  ) 


APPENDIX  C 


PROGRAM  LISTING 


This  appendix  contains  the  c-omplefce  program  listing. 
All  subroutines  are  stored  on  disc  6.  For  execution,  it  is  only 
necei sary  to  compile  the  MAIN  program.  The  tvo  subroutines  M0CP 
and  M0CB  are  stored  in  the  LOCAL  mode. 
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MAIN  PROGRAM 
COMMON  STATEMENTS 

COMMON  NU,  CONST, THETW, PFLRT,  PI , X3A, Y3A,E3A, V3A,T3A, P3A,R3A,  VSLP3, 
ITPTA  RPHA*7FPTA  PSITA 

COMMON  GAMMA,GC^RGAS,DI AM, RHOP,  PR,  CPART 

10  CONTINUE 
PAUSE 

CALL  DATSW  (10,K1Q> 

GO  TO  (20,3ft),  K10 
20  CALL  FLO V/P 

no  to  io 

30  CALL  EXIT 


o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  n  o  o  r>  o  o  o  .">  r>  ri  ;■>  o  n  r>  o  o  o  o  o  o  o  o  o  o  .o 
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SUBROUTINE  FLOWP 


THIS  PROGRAM  CONTROLS  THE  LOGIC  NECESSARY  TO  CALCULATE 
THE  FLO'/  FIELD  IN  A  SUPERSONIC  AXI  SYM.METRI  C  (NU«1)  OR 
TWO-DIMENSIONAL  (Nti-0)  NOZZLE  BY  THE  METHOD  OF  CHARACTER  IS 

THE  ARRAY  SLV  IS  USED  TO  STORE  THE  I  NIT. AL- VALUE  LI  HE.  . 

THE  FIRST  SUBSCRIPT  IS  THE  POINT  NUMBER,  THE  SECOND  DENOTE 
THE  Fl'ID  PROPERTY  DEFINED  BY  THE  FOLLOWING 

1  ;;  position 

2  Y  POSITION 

3  GAS  MACH  NUMBER 

4  GAS  VELOCITY 

5  GAS  FLOW  ANGLE 

6  GAS  P.RESSURF 

7  GAS  DENSITY 

3  PARTICLE  VELOCITY 

9  PARTICLE  FLOW  ANGLE 

10  PARTICLE  DENSITY 

11  PARTICLE  TEMPERATURE 

12  FARTICLE  STREAM  FUNCTION 

THE  SAME  SCHEME  IS  USED  FOR  THE  TWO  CHARACTERISTIC  ARRAYS 
CHAR A  AND  CHARL.  ,  TWO  ARRAYS  ARE  NEEDED  SINCE  IT  IS  NECESS 
TO  SAVE  ONE  LEFT  CHARACTERISTIC  WHILE  THE  NEXT  ONE  IS  BE  IN 
CALCULATED. 

**  M0Tc:  ** 

*i  HE  INPUT  PARAMETER  iJPTSL  GOVERNS  THE  SIZE  OF  THE  THREE 
ARRAYS,  THESE  MUST  CE  5LV(NPTSL, 12 ),  CHARA(2*NPYSL, 12) 

AND  CHARL(2*HPT5L/12). 


TYPE,  DIMENSION  AND  COMMON  STATEMENTS 

n  jr a  l 

DIMENSION  S LV(5,  12 ),  CHARAdS,  12),  CHARL (10,  12) 

COMMON  N*l,  CONST,  THETW,  PFLRT,  PI  ,  X3A,  Y3A,  E3A,  V3A,T3A,  P3A,  R3A,  VF  '3, 
1TP3A, RP3A, TFP3A, PS  I 3A 
COMMON  GAMMA, GC, RGAS . 0  I  AM, RHOP, PR, CPART 

SET  PROGRAM  CONSTANTS  A NO  IMJATILIZE  COUNTERS 

GC=0.S,1 
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n  n  c >  no 


-  A3-5  - 


C 

40  J»!+l 

GO  TO  GO 
C 

50  J^MPTSL+HPTSL-l 
ti=j 


THE  CC'JHTER  r.U  RETET1IUES  WETHER  the  left  characteristic 
•'.HI  HG  CALCULATEO  STARTS  FRO  1  THE  INITIAL-VALUE  LIME  OR  FRO 
THE  CEHTERLIHE 


5R  i:*>i+i-;;?tjl 


0 


r 

> 


r* 


LT)p  FOR  CALCULATING  THE  PROPERTIES  A  LOUR  A  LEFT  CHARACTER 


!'■ 

80 

0 


10C 


n°  2nn  L»1,J 
l;  =  L+l 
•*=L+R 

:.=l-i 

!  F  (f! }  V' ,  7f',  12rt- 
contipt 

IF(KK)  RR,  8n,llR 

IN0EX=2 

JP=I+1 

no  100  10=1,12 
CHARA(L,IP)=SLV(JP,!P) 
CONTINUE 
no  TO  120 


CALCULATE  FLOW  PROPERTIES  T’lE  C“  JTERLI  ME 


110  ICL=2 

vf.lpi=c:;a  >.1(2,8) 

VELP2=r/:/.  IL(1,  C) 

TP1-C:!A.RL(2,3) 

TP2  =  CHARL  ( 1,  0 ) 

.RP1®C!IARL(2, 10) 
rp2=o;:arl(i,  io) 

TEP1*CMA  'L(2, 11) 

TEP2=R:;AR  L(  1,  11) 

PS  1 1-CHARLC2, 32 ) 

:>3 1  2=c:A">L(  1  12) 

CALL  •  •''CO  ( 2 ,  CHAR L < 2 ,  1 ) ,  PH ARL  (1,  1 ) ,  Cl !AR L ( 2 , 2 ) ,  CHAR L ( 1,  2 ) ,  CHARL (2,3) 

1,  CHARLC 1,  3),  CHAP. L (2, 4),  CIARL(  1,  4 ) ,  CMARL( /,  5),  CHARLd,  5),  CHARLd,  6) 

2 ,  CHAR L(  1,  C ) ,  C! !AR L ( 2,  7 ) ,  C! I AR L O  ,  7 ) ,  VE L PI,  VE LP2,  TP1,  TP 2,  R  PI,  R P2, 
;>TEP1,TEP2,PS!1,PS!2) 


A3~6 


CO  TO  If,?'. 


12"  co:jti:"f 

CO  TO  (130, UD,  I  CL 
I3r  i>L 
::=L 

CO  Tn  15*' 


\« 


c 


r 


14° 

15T  VFLPl'O'LV'LO'.,^) 
VELP2  =  .riA'.A(f,f) 
TPl=Cl'A,'L{!»#  0) 
TP2=CHATACI,0) 

h»»i»c::atl(‘'.,ia) 

CP2=CMATA(:;,  ID 

TEP1=CMATL(K, 11) 
TEP2=C;!A.TA(U/  11) 
?3!i»C:!.V'.L(!'.,12) 

PS  I  2=CHA  TAO!,12> 

\*fsi  L— j 

IF (:;,:)  ISC,  170,  17C 


CALCULATE  FLJ!;  PTOPETTIFS  AT  .VI  li.TETIOT 


POI 


r.A 

LL 

13  CP 

1), 

ci:.\ 

TA(H, 

2), 

r*'i  a  n  *  f  m 

V/  l/V  v/ll  . 

3TEP1, 

tep:>. 

rr5 

TO 

in" 

( c;  I  AT  L(K,  U,  Ci’ATA  ( IS,  1 )  ,  CM  AT  L  ( K,  2 )  ,  CHAT  A  ( iJ,  2  )  ,  C!  I  AT  L  ( X  ,  3 
3),  CHATL( ;;  ,4),  CHAT  A  (  U,  4 ) ,  Ci  I  AT  L  (  K,  5 ) ,  Cl  I  AT  A  (  H,  5 )  ,  CHAT  L  (  K,  R 
3  )  ,  CH  AT  L  ( K  ,  7  )  ,  C! !  AC!  A  ( 1 1,  7),  VELP1,  VCLP2,TP1,TP2,TP1,TP2, 
PSP.,  PS!  2) 


C 


r 

\# 

c 


r 


CALCULATE  f LT l  PTOPETTIFS  OH  THE  UOIZLF  ’..’ALL 


17"  CALL  ’.0  CO  ( 1,  C’ !  AT  L  (!'.,!),  CHAT  A  ( !J,  1 ) ,  C"  AT  L  (  K  ,  ? ) ,  CH  AT  A  ( !J,  ? ) ,  CHAP.  L  (  K  ,  3 
D  ,  C**  AT  A  ( 3  )  ,  CH  ATL  <  K,  4  )  ,  CHATA  ( M, 4 ) ,  C!  'AT  L( K, 5 ),  CHAT  A  (M,  5),  C!!ATL(K,  P 
2  ),  C,!ATA( G),C!!A?L{T,7),  CMATAOJ,  7), VE  LP1, VELP2,  TP1,  TP2,T PI,  RP2, 
3TEPl,TEP2,PSn,PSI  2) 


ST-iTE  THE  SOUJTIOU  0TTAIi.,rO  FTO.H  SUTROUTbHE  HOC  !  0 
THE  AUX ILL! ARY  ATTAY 


IS''  CHATACH,  1)  =  X3A 
CHATAC1,  2)=Y3A 
CHATA  Cl,  3  )  =  E3A 
CHATA  (■’*.,  4  )=V3A 
CHATA  C  I,  5)=T3A 


A3-7 


CHA1AC1  /6)*P34 
CHARA(M/7)=R3A 
CHARA(M/3)«VELP3 
CHARA(M/3)=TP3A 
CHAR  ACM*. 10)=RP3A 
CHARACM/ 11)“TEP3A 
C!IARAC1/12)*PS!3A 
200  CONTI  HUE 
210  CONTINUE 


MOVE  OATA  FROM  AUXILLIARY  ARRAY  TO  THE  LEFT 
CHARACTERISTIC  ARRAY 


DO  24C  K-l/MM 
no  22n  11=1/ 12 
CHARLCX/ LL)=CHARA(K/ LL) 


22C  CONTI  nun 

I F C I -I OUT)  240/230/230 


PRINT  OUT  SOLUTION  ALONO  7HF  LEFT  CHARACTERISTIC 


230  CONTINUE 

I F (L I NE-62)  234/232/232 
232  WRITE  (1/1008) 

WRITE  (1/1006) 

Ll NE=9 

234  XCM=C!IARL(K/1)*METER 
YCM»CHARLC</2)*METER 
TO  EG=  CHAR  L ( K / 5 ) *R  AO 
RKG=CIIARL(K,  7)/GC 
PRES=CMARL(K/S)/(GC*.METER*FiETER) 
TE'1P=CHARL(fC/6)/(RGAS*CHARL(K/7)) 

TPnEG=CHARL(K/9)*RA0 

WRITE  (1/1004)  XCM/YCM/CHARL(K/3)/CHARL(K,4)/TnEG/PRES/RK! 
1/TEMP 

WRITE  (1/ 1007)  CHARL(K/8), TPOEG/ CHARL(K/ 10)/ CHARL(K/ 11) 

LI ME=L I NE+2 
240  CONTINUE 

IF(I-IOUT)  260/250/250 
250  WRITE  (1/1005) 

LINE=LliIE+3 

STORE  WALL  POINT  IN  THE  LEFT  CHARACTERISTIC  ARRAY 


260  im-M.M+1 

HO  2  70  LL=1/12 
charl(m:.*/  ll)=ciara(.m.M/  ll) 


on 


-  A3-8  - 


2  7n  COUTI 

I F  < I -I  I  AX)  30,30,230 

2  30  RET'JRi! 

FORMAT  STATEMENTS 

10C1  FORMAT  (4FIO.rt) 

3,^02  FORMAT  (  I  2, 3X,  I  2, 3X,  I  2,  SX,  1 1) 

1003  FORMAT  ( 7F1H. C/5F1Q.  0) 

lrt04  FORMAT  (3(F  10. 5,  IX),  2(F  10, 4,  IX) ,  Fll.  4.,  3X,  F1P.  4,  IX,  F 10. 2,  8X,  1  GAS’  ) 
1^5  FORMAT  <//) 

iron  FORMAT  (2GX,  ‘MACH’,  2AX,  ‘ F LO’./’/GX,  ' X' ,  10X,  'Y*,  7X,  *  NUMBER’,  4X,  ‘VELOC 
IITY’^X,  'ANGLE ',  ;x,  'PRESSURE',  5X,  '  DEflSI  TY' ,  3X,  ‘TEMPERATURE  * //4X, 

2  '  (CM) ',  7X,  '  (CM)  ',  17X,  1  (M/S EC ) '  ,  6X,  '  (OEG)  1 , 4X,  '  (KG/SQ  CM)  ',  3X, 

3  *  CKR/CTJ  M)*,4X,  '  (DEG  K)  »//) 

3. nri7  FORMAT  ( 33X,  2(F  Iff.  4.,  IX),  14X,  FIR. 4,  IX,  F10. 2,  SX,  *  PARTI  CLE  * ) 

10R3  FORMAT  (////) 

C 

f.:io 


onoooooo 


1 


*•*  A3~9  ** 


SUB  ROUT  I  flE  COEF  P  ( Y,  E,  V,  7,  P,  RHO,  VELP,  TEP,  RHOPP,  TP,  S,  H,  Q,  APU,  APV, 
IARU,  ARV,  CE,  G,F,  FP,  GP,BP,  C,  OP) 


THIS  SUBROUTINE  CALCULATES  THE  COEFFICIENTS  AT  EACH  POINT 
UMICH  ARE  USED  IN  SOLVING  THE  DIFFERENTIAL  EQUATIONS  FOR  THE 
METHOD  OF  CHARACTERISTICS  SOLUTION  TECHNIQUE 


TYPE,  DIMENSION  .‘40  COMMON  STATEMENTS 

REAL  K GAS, MU, MU1 
C 

COMMON  IIU,  CONST,  THETW,  FF  LRT,  PI  ,  X3A,  Y3A,  E3 A,  V3A,  T3A,  P3A,  R3A,  VEI.P3, 
1TP3A,RP3A,TEP3A, PS  I 3A 
COMMON  GAMMA, GO, R,D l  AM, RHO P, PR, CPART 
C 

SOS*V/E 

n*i.a/so.RTCE**2-i.a) 

A*ATAN(D) 

STMA*S I N(T-A) 

STPA«SIM(T+A) 

CTMA*COS(T-A) 

CTPA*COS(T+A) 

S*STMA/CTMA 

H*STPA/CTPA 

Q  roS(A)/(SIN(A)*RHO/GO*V**2) 

UP»VELP*COS(TP) 

VP“VELP*SIN(TP) 

DU*V*COS(T)-UP 
DV=*V*SI  N(T)-VP 
TEMP*P/(R*RHO) 

MU»0.n,0000143*SQRT(TEMP)/(l.a+105.a/TEMP) 

KGAS=10no,fj*MU 

REY*RHO*OI AM*A8S( V-VELP)/ (MU*GO) 

CX*0.U.8+28.P,/(REY**0.a5) 

NU1=2,P.+  P,6*SQRT(RLY)*PR**0,3.33 
CP»12.P*KGAS*NU1/(MU*CX*REY) 

R=  ( GAMMA-1 .  n.)*(2.n*CP*CTE P-1  EM P ) / 3 , OU*DU+OV*OV ) 
AP=G.7.5*CX*RHO*ABSCV-VELP)/(D!  AM*RHGP*GO) 

APU=*AP*OU/UP 

APV=AP*OV/UP 

ARU»AP*RHOPP*DU 

ARV=AP*RHOPP*DV 

CE*  2  .  P,*A  P*  CP*  (TEMP-TE  P )  /  ( 3 . 0.*U P* CPART ) 
F*-AP*GO*RHOPP/(RMO*V*E*STPA)*  ( 1.  O.+  B/  ( SOS* SOS)  ) 
G»-AP*On*RHOPP/(RHO*V*E*STnA) * (1,0+0/ (SOS*SOS ) ) 

IFCY-Q.noonooi)  20,_20,10 

10  STY*SIH(T)/Y 
F-F+STY/(E*STPA> 

G«G+STY/(E*STMA) 

20  FP*-AP*RHOPP*GO*VELP/(V*V*RHO)*(COS(TP)*H-SJN(TP)) 
GP«-AP*RHOPP*GO*VELP/<V*V*RHO)*(COS(TP)*S-S!N(TP>) 

BP«RHO*V/GQ 

C“SOS*SOS/GO 

DP«AP*8*RH0PP/(C*V*C0S(T) ) 

RETURN 

END 


SUBROUTINE  HOCP  (XI, X2, Yi, Y2,£l, E2, VI, V2, Tl, T2, PI, P2,R1,R2> 
lV£tP2,V£LP2*TPl,TP2,RPl,RP2,TEPl,TEP2,PSIl,PS!2) 


THIS  SUBROUTINE  IS  A  METHOD  OF  CHARACTERISTICS  SUBPROGRAM 
USED  TO  CALCULATE  THE  FLO’./  PROPERTIES  AT  AN  INTERIOR 
MESH  POINT.  . 

THE  GASDYNAMIC  MODEL  IS  ROTATIONAL  AND  A  PERFECT  GAS  IS 
ASSUMED,  ,  THIS  SUBROUTINE  USES  THE  MOOI F I ED-EULER 
PREDICTOR-CORRECTOR  TECHNIQUE,  WHERE  AVERAGED 
COEFFICIENTS  ARE  USED.  . 

POINTS  1  AND  3  ARE  ON  THE  RIGHT  CHARACTERISTIC,  POINTS 
2  AND  3  ARE  ON  THE  LEFT  CHARACTERISTIC,  POINTS  5  AND  3  ARE 
ON  THE  GAS  STREAMLINE,  POINTS  4  AND  3  ARE  ON  THE  PARTICLE 
STREAMLINE  WHERE  POINT  3  I S  THE  DESIRED  SOLUTION  POINT.  - 


TYPE,  DIMENSION  AND  COMMON  STATEMENTS 

COMMON  NU, CONST, THETW, PFLRT, PI , X3A, Y3A, E3A, V3A, T3A, P3A, R3A, VELP3, 
1  TP 3 A,  R  P3 A ,  TE  P3A ,  PS  I  3 A 
COMMON  G, GO, R, D I  AM, RHOP, PR, CPART 
l*-l 
ICL-1 


CALCULATE  THE  COEFFICIENTS  AT  POINT  1 

CALL  COEr  m,  El,  VI,  Tl,  PI,  Rl,  VELP1,  TEP1, RP1,  TP1,  SI, H,  Ql,  API),  APV, 
1ARU,  ARV,  v  G1,F,  FP, GP1,B,  C,  D) 

AVERAGE  THE  COEFFICIENTS  FOR  THE  PREDICTOR 


S13  B1 
Q13-Q1 
G13-G 1 
GP13-GP1 


CALCULATE  THE  COEFFICIENTS  AT  POINT  2 

CALL  COEfP  <Y2,E2,V2,T2,P2,R2,VELP2,TEP2,RP2,TP2,S,H2,Q7,APU,APV, 
1ARU, ARV, CE, G2, F 2  JP2,GP,B,C,D) 

I  F  (Y2-0. 0,000001)  70,70,80 


CENTERLINE  APPROXIMATION 


ooooo  o  ooo  noon  non 


A3 -11  - 


70  I CL-2 


AVERAGE  COEFFICIENTS  FOR  THE  PREDICTOR 

1 

80  H23*H?  | 

Q2  5*Q2  I 

F23-F2  1 

FP23*FP2  | 

-------  j 

SO  1*1+1  | 

INTERIOR  MESH  POINT  SOLUTION 

X3A* ( Y2-Y1-X2*H23+X1*S13)/ (S13-H23) 

Y3A*Y1+S13* (X3A-X1) 

GC  TO  (159/140),,  I  CL 

SOLUTION  WHEN  POINT  2  IS  ON  THE  CENTERLINE 

140  P3A^{-2.Q*U+Q23*P2+2.n*Q13*Pl-2.P.*G13*(Y3A-Yl)-2.0+GP13*{X3A-Xl> 
l-F23*{Y3A-Y2)+FP25*(X3A-X2))/(  2. 0*0.15+023) 

T3A*T1+013*  ( P3A-PD+G13*  (  Y3A-YD+GF13*  ( X3A-X1)  ! 

F 2*F  2+T3A/Y3A  I 

I CL*1  ! 

GO  TO  210  I 

150  P3A*(Q23*P2+013*P1+T2~T1-F23*£Y3A-Y2)+FP23*CX3A-X2)-G13*(Y3A-Y1)~  [ 

inP15*(X3A-Xl) )/(023+Q13) 

T3A-T2-Q2  3* ( P3A-P? ) -F2  3* ( Y3A- Y2 ) +F  P23* ( X3A- X? )  | 


CALCULATE  THE  POSITION  OF  POINT  5  FROM  THE  GAS  i 

STREAMLINE  CHARACTERISTIC  EQUATION.  ,  { 

210  SXY-1 

JX12«X1-X2 

IFCABS (0X12) -0.0.091)  230/240,240 
230  DX0Y=(X2-X1)/(Y2-Y1) 

DX12*Y1-Y2 
I  XY-2 
GO  TO  250 

240  0YDX»(Y1-Y2)/DX12 
250  CONTINUE 

I  F  ( I  )  2  60,260/270 
260  T5*T3A 

270  TT35* (SI  N(T3A)/COS(T3A)  +S I  N(T5 )/C0S (T5)  )*0.  5, 


i 


i 


s- 


280 


290 


GO  TO  (280,290),  l  XY 

X5« (Y3A-Y2+  X2*DYDX-X3A*TT35 ) / (DYDX-TT35) 
Y5*Y3A+TT35* (X5-X3A) 

OX-X5-X2 
GO  TO  300 

Y5*  (Y3A+TT35*  (X1-X3A)-TT35*DXDY*Y1 )  /  ( 1. 0,-TT35*DXDY) 

X5»X1+(Y5-Y1)*DX0Y 

OX-Y5-Y2 


C 

C 

C 

C 


CALCULATE  THE  DERIVATIVES  OF 
POINTS  1  AND  2 »  . 


THE  PROPERTIES  BETWECM 


C 

C 

C 


300  SM12«(E1-E2}/0X12 
SV12* (V1~V2)/DX12 
ST12* (T1-T2) /DX12 
SP12*(P1-P?)/0X12 
SR12*(R1-R2)/DX12 
SPV12* (VELP1-VELP2 ) /DX12 
STP12® (TP1-TP2 ) /OX12 
SRP12»(RP1-RP2)/DX12 
STE12® (TEP1-TEP2)/0X12 


C 

C 

C 

C 


502 

304 


306 


LINEARLY  INTERPOLATE  FOR  THE  PROPERTIES  AT  POINT  5 


E5*E2+SM12*DX 

V3*V2+SV12*OX 

T5*T2+ST12*DX 

P5-P2+SP12*DX 

R5»R2+SR12*DX 

VELP5aVELP2+SPV12*DX 

TP5-TP2+STP12*DX 

RP5=»RP2+SRP12*DX 

TEP5aTEP2+STE12*OX 

I F ( I )  302,302,30b 

CALCULATE  THE  POSITION  OF  POINT  4  FROM  THE 
STREAMLINE  CHARACTERISTIC  EQUATION.  . 

TP3A=TP5 

TP4=*TP5 

TT34=( SI  N(TP3A)/COS (TP3A)+S I H(T P4)/COS  (TP4))*0.5, 
GO  TO  (306,307),  I  XY 

X43 ( Y3A-Y2+X2*OYOX-X3A*TT34)/ (DYOX-TT34 ) 
Y4*Y3A+TT34* ( X4-X3A ) 

DXX*X4-X2 


PARTI CLE 


O  O  O  O  Ci  O  O  O  O  o  o  o 


c 

307  YWY3A+TT34*(X1-X3A)-TT34*DXDY*Y1)/(1.0-TT34*RX0Y) 
X4*XI+CY4-Y1)*DXDY 

DXX*Y4-Y2 

C 

C  LINEARLY  INTERPOLATE  FOR  THE  PROPERTIES  AT  POINT  4 

£ 

308  E4»E2+SM12*OXX 
V4»V2+SV12*OXX 
T4»T2+STI2*OXX 
P4«P2*SP12*OXX 
R4»R2+SR12*DXX 
VELP4»VELP2+SPV12*DXX 
TP4»TP2+STP12*DXX 
RP4»RP2+SRP12*DXX 
TEP4»TEP2+STE12*DXX 

PS  1 4*PS  I  2+PI  /PFLRT*  ( 0. 5.*  (RP4*  VELP4*COS  (TP4 )  «-RP2*VELP2*COS(’rP2)  ) 
l*(Y4**2-Y2**2)~(Y4*Rr'  VELP4*SI N{TP4)+Y2*RP2*VELP2*SI N(TP2) ,*(X4 
2-X2) ) 


CALCULATE  THE  COEFFICIENTS  AT  POINT  5 


CALL  OQEFP  (Y5,E5,V5,75,P5,R5,VELPS,TEP5,RP5,TP5, &PU, APV, 
1ARU5,  ARV5,  CE,  G2,F,FP,  GP,B5,C5,05) 

CALCULATE  THE  COEFFICIENTS  AT  POINT  4 

CALL  COEFP  (Y4,E4,V4,74,  P4,R4,VELP4,TEP4,RP4,TP4,S,H,Q,APU4,AP74, 
1ARU,ARV,CE4,G2,F,FP,GR,B,C,D) 

IFU)  >2fi,320#33« 

AVERAGE  COEFFfCh  TS  FOR  THE  PREOICTOR 

32P  C35*C5 
335=B5 
035*05 
ARU55“ARU5 
ARV35*ARV5 
APU43*APU4 
APV43-APV4 
CE43-CE4 
GO  TO  54 U 

AVERAGE  COEFFICIENTS  FOR  THE  CORRECTOR 
330  S 35*{ B3*BS )*P« 5, 


-  - 


C 

c 

c 


035*0. 5*(03+f)5) 

ARU35  =  0«5,*  (ARU3+ARU5) 

ARV35=Q.  5,*  (ARV3+ARV5) 
APU43=0, 5> (APU4+APU3) 
APV43=0.5* (APV4+APV3) 
CE43  =  0.5,*(CE4+CE3) 


r  LCULATF.  THE  VELOCITY  AND  DENSITY  FROM  THE  STREAMLINE 
^.'MPATIBI  LITY  EQUATIONS 


340  V3A*V5+  ( P5-P3A-ARU35*  ( X3A-<^5)  -ARV35* ( Y3A-Y5 )  )/B35 
R3A=R5-s-(P3A-P5)/C35-D35*(X3A-X5) 

TEMP*  P3A/ (R*R3A) 

S OS* CONST* SQRT( TEMP) 

E5A*V3A/SOS 

UP3»VELP4*COS(TP4)+APU43* (X3A-X4) 

VP3*VELP4-  SI NCTP4)+APV43* (X3A-X4) 

VELP3»SQRT(UP3**2+VP3**2) 

YDOT3-VP3/UP3 

TP3A-ATANCYD0T3) 

TEP3A=TEP4-*-0E43*  (X3A-X4) 

PS  I 3A=PSI 4 

RP3A*  ( PFLRT* ( PS  I  3 A- PS  1  2 ) / PI  -RP2*VELP2*COS  (TP2  )*Q.  5,*  ( Y3A**2-Y2**2) 
1+(X3A-X2)*Y2*RP2*VELP2*S!NCTP2))/<0.5>'JP5*(Y3A**2-Y2**2)-Y3A*VELP3 
2*SI H(TP3A)* (X3A-X2) ) 

I F ( I }  350, 350,410 
C 

C  CALCULATE  THE  COEFFICIENTS  AT  POINT  3 

C 


350  CALL  COEFP  (Y3,?  ,  E3A,  V3A,  T3rt,  P3A,  R3A,  VELP3,  TEP3A,RP3A, TP3A,  S3,  H3, 
1Q3, APU3, APV3, ARU3, ARV3, CF3, G3, F3, F  P3,GP3,B5,C3,D3) 


C  AVERAGE  COEFFICIENTS  FOR  THE  CORRECTOR 

C 

H23*  (M2  *H3)*0, 5. 

Q23*(Q2  +  Q3)*0.5, 

F23*  (F2+F3)*0.  5, 

$1$»(S1+S3)*0.5. 

Q13* (Q1+Q3 )  *0, 5, 

013'  (Gl-r,3)*0,5. 

FP23=0,5*<rP2+FP3) 

GP13=0.5* (GP1+GP3) 

GO  TO  90 
C 

410  RETURN 
END 


i 


_  a-s-IC  _ 


C 

C 


crvr.ori  ::u,  const,  thetu,  pflrt,  Pi, X3A, Y3A,C3A,v3A,T3A,  p 3A,R3A, vel?3. 


1TP3A,PP3A,TEP3A,  PS  I  3 A 
CO;  7 10  M 
l=-l 
r.o  to 


c 

c 

c 


o,  o'1 ,  a,  o  f  a;  t,  an  op,  pr,  c  part 
( lr, 20 ),  10 

REDEFINE  DATA  TRANSFERRER 
-./ALL  PO I  :JT  SOLUTION 


10  X 5» XI 
Y5»Y1 
E5a£l 
V5=V1 
T5=T1 
P  S33  PI 

r>5  =  0  1 

VELP5=VELP1 
TP5=TP1 
RP5=RP1 
TF.?5=TEP1 
PS  I  5^  PS!  1 

TAMTU=$ ! !!  ( T! iF.T. J)  /CO S  C THET’ ') 

X4=X1 
Y  4=Y1 
L4  =  n 
V4=V1 
T4=T1 
P4=P1 
R  4  =R  1 

VFLP4=VELP1 

TP4=TP1 

RP4=RP1 

TEP4=T-E?1 

PSI4-PSI3 

GO  TO  4  0 

..F.RLFINE  data  transferred 
centerlimf  soltjt s o:i 


THROUGH  CALL  STATE" 


FOP. 


THROUGH  CALL  STATEMENT  FOR 


2° 


X5=X2 
Y5=Y2 
E5=i:2 
V  5=V2 
75rT2 
F5'»p2 


o  c:  c 


/3  -l6 


0  5=R2 

VELP5=VELP2 

TP5=TP? 

•i  {>  5  =  o  r>2 

fEP5=TEP2 
PS  !  5=  PS  I  2 
X4=X2 
Y4=Y2 
r4  =  E2 
V4=V2 
T  4  =  T2 
P4=P2 
'’.4=°.  2 

7ELP4=VELP2 
TP4=TP2 
P  P  4 = 0  P2 
TFP4=TEP2 
PS  I 4  =  PS 1 2 

CALL  C^FPP  (Yl/El,Vl,Tl,Pl,Rl,VELPl/TEPl/RPX,TPl,Si,H,0.1,APU,APV, 
1 ACJ,  AR7,  CE,01,r,F P,  GP1,  0, C, D ) 

SI  3= SI 
:'.l3=Qi 
fi 13=01 
C.  PI  3=0  PI 

fin  TO  IS 

4  G  CALL  C.OEFP  (Y2,  E2  ,  V2,  T2,  P2,  R2«  VELP2,T£P2,RP2,  TP2,  S,  HI,  Q2,  APIJ,  APV, 
1ARU,  A HV,  CF,  32/F2/FP2/0P/  0,0,0) 

M23=H2 
Q23=Q2 
F  2  3=  F  2 
F  P23=i:  P2 
DC  1=1+1 

00  TO  (IP*, ISO),  IP. 

’.’ALL  POINT  SOL'JTIn*l 

IOC  X3A=  ( Y5-Y2+I!23*X2-X5*T»V!T’ «)/  (H23-TA*'T'.J) 

Y3A=Y?+!!23*  (X3A-X2) 

t3A=T':ft\: 

P3A=  (02  P2+T2-T3A-F23*  ( Y3A-Y2 )  +FP23*  CX3A-X? )  )/')2  3 

00  TO  2  lr' 


OF  !Trr  _l  ME  SOLMTIOIi 


ISC  X5A=X1-Y1/S 13 
Y34=P.n 


*  A3«l?  • 


T3A=r.r-, 

P3A=(ni3*Pl-Ti+r.l3*Yl-rjP13*(X3A~Xi)  )/ni3 
210  CONTINUE 

I  F  ( I  )  31P, 31", 330 

3 1°  CALL  COEFP  <Y5,E5,V5,T5,P5,R5,VELP5,TEP5,RP5,TP5,S,H,Q,APU,APV, 
J.AOU5,  A0V5,  CE,  fi2,  F,  ~  P,  C-P,  3  5,  C5, 05 ) 

CALCULATE  T NO  COEFFICIENTS  AT  POINT  4 

CALL  CHEF P  ( Y4,  E4, V4,  T4,  P4, R4,  VELP4,  TEP4,  RP4,  Tp4,  S, !!,  Q,  AP'JU,  APV4, 
lAI'J,  AO  7,  CO  4,  52, F, F P,  CP,  3,  C,  0) 

C35=C5 

3  35=B5 

D35*05 

ARU35-AOU5 

ARv/35-ARV5 

AP’J43=APU4 

AP743-APV4 

CF.43=CE4 

SO  T  34 P 

3 5C  P.3S=(?.3+C5)*n.5 
035=*  C3  +  C5W.5 
035=*. 5* (03+05) 

AOU 3 5=  0.5*  (A  VJ 3 + ARU 5 ) 

; \V35= 0. 5*  </ OV3+ARV5) 

A?U43=r’»  5*  (APM4+APU3) 
t D743=  F. 5* ( APV4+A  PV3) 

CE45=r.5*(cr4+CF3) 

3  4  0  7  3  A = V  5  +  ( P  5  -  P  3  A  -  -V  VJ  3  5  *  <  X  3  A  -  X  5  )  ~  A  0  V  3  5  *  ( Y  3  A  -  Y  5  >  )  /  R  3  5 
?3A=r.  5  +  ( r*3A-P5  )/C35-035*  (X3A-X5 ) 

TF'IP=P3A/  (0*0 3 A) 
s^G=co;;sT*scrf(TE’iP) 

F3A=73A/50S 

UP3=VELP4*C0S(TP4)+AP'J43*  (X3A-X4) 

VP3=VELP4*S I  :(TP4)+APV45*(X3A-X4) 

VELP3=SQ  1T(UP3**2+y?3**2 ) 

YOOT3=VP3/UP3 
T  n  3  A = A  T  A!  I  ( Y  00  T3 ) 

TEP3A=TE°4+CF43* (X3A-X4) 


PS  I 3A= OS  1 4 

342  0?3A=  (PFLOT*(Pr/l  3A-PSI  2 )  /  PI  -0P2*7ELP2*COS  (TP2 )  *Q.  5*  ( Y!A**2-Y2**2) 
1+ (X3A-X2)*Y2*RP2*VELP2*S i N(TP2) )/ ( 0. 5*UP3* ( Y3A**2-Y2**  2) -Y3A*VELP3 
2*5  !■  ii(TP3A)*  (X3A-X2)  ) 

SO  TO  34C 

344  0  P3A~  { pf:  LOT*  PS  1 1/PI -OP1*V£LP1*COS (TP1) *0. 5*Y1**2+ (X1-X3A) *Y1  i* 


o  o  r> 


3LVELP1*SI  ::<TP1> )/  (0.5*UP3*Y1**2) 

346  CONTI  MU': 

!F(!)  35^, 350, 410 

35r  CALL  COEFP  (Y3A,  E3A,  V3A,  T3A,  P3A,  U3A,  VELP3,  TEP3A/RP3A/TPSA/  S3#H3/ 
1Q3,  APIJ3,  APV3,  AR03, A0V3,  CF.3/  G3,  F3,  FP3,  GP3,  S3,  C3,  D3  ) 

CO  TO  (300,370),  tO 

CENTER  LIME  A  PPM  OX  t '  tAT !  0  N 

370  STY=SIM(T1)/Y1 

0=1,  0,/S  OPT  ( E3A**2-1.  ”) 
a=--ata::(o) 

STMA=S!N(T3A-A) 
n  3=0  3+ STY / ( E 3A+STMA ) 

00  TO  4-no 

300  •!23=(!!2+!'3)*C,5 
0.27/*  (02  +  03 )  *0,  5. 

F23*  F 2+P3 )*Q,  5 
FP23=0.5*(FP2+FP3) 

00  TO  (0C,4np),  |  0 

4 OP  S 13=(ol+33)*0. 5 
n.13=  (Ql+03  )*  0, 5 
013«<G1+G3)*C.5. 
npi3=o05*(npi+np3) 
no  to  90 
410  HE  Tin  I! 

HMD 


t 


-  A3-19  - 


C  SAMPLE  DATA  CASE 


1.4, 

29.2,7 

0.7, 

0.0,005 

2000,0. 

500. a 

0.8,972 

3 

35 

0 

1 

187.9,731 

0.0, 

0.005485 

1.0,736 

330.3.1 

15,0. 

1296100. a 

123.7,6 

15.0, 

76.7,112 

2  98.0.5 

l.ft 

18  7.9,731 

0.0,0002 

0,0.02  7425 

1.0,736 

330.3,1 

7.5. 

1296100. a 

123.7-6 

7.5. 

76.7,112  8 

298.0.5 

0.25 

187.9,731 

0.R0004 

0,0. 

l.a?36 

330.31 

6.0, 

1296100.0, 

123.7,6 

o.n. 

76.7,112 

298.0.5 

0.0. 

' 

